#include "TH1D.h"
#include "TH1D.h"
#include <vector>
#include <iostream>
#include "TROOT.h"
#include "TMath.h"
#include "TCanvas.h"
#include "TRandom.h"
#include "TF1.h"
#include "TFile.h"
#include <TStopwatch.h>

void main(int random_seed) {

	const double Pi = TMath::Pi();
	const double pi = Pi;
	TRandom *myRndm = new TRandom(random_seed);
	double Psi[4] = {0,0,0,0};
	double v2_input=0.1, v2_measured=0;
	const int N_Particle = 500;


	for (int ih = 1; ih <4; ih++){
	Psi[ih] = (myRndm->Uniform(-1 * pi / ih ,  pi / ih )) ;
	}


                                  // coeffi   const   1th order           2th order           3th order // 
	TF1 *dNdphi = new TF1("dNdphi", "[0] *(   [1]+  2*[2]*cos(1*(x-[3])) +2* [4]*cos(2*(x-[5])) +2* [6]*cos(3*(x-[7])))", 0,2* pi);
	dNdphi->SetParameter(0, 100); 	 //coeffi
	dNdphi->SetParameter(1, 1);  	  // const term
	dNdphi->SetParameter(2, 0); 	 // directed flow magnitude  (v1)
	dNdphi->SetParameter(3, Psi[1]); 	 // directed flow direction  (psi 1)
	dNdphi->SetParameter(4, v2_input); 	 // elliptic flow magnitude (v2) 
	dNdphi->SetParameter(5, Psi[2]); 	 // elliptic flow direction (psi 2) 
	dNdphi->SetParameter(6, 0); 	 // tri flow magnitude (v3)
	dNdphi->SetParameter(7, Psi[3]); 	 // tri flow direction (psi 3)

	double phi[N_Particle];
	for (int i=0; i<N_Particle; i++){
		phi[i] = dNdphi->GetRandom(); 
	}	
// generation done


// find v2 magnitude and direction
				// cn  1th 2th  3th
	double An[4] =  {0, 0,  0, 0};
	double Bn[4] =  {0, 0,  0, 0};

	int ih=2; // for (ih = 0; ih < 4; ih++);

	for (int ip=0; ip<N_Particle; ip++){
		An[ih] += TMath::Cos( ih* phi[ip] );
		Bn[ih] += TMath::Sin( ih* phi[ip] );
	}
	
    double C_Psi[4] = {0,0,0,0};
	C_Psi[ih] = TMath::ATan2(Bn[ih],An[ih]) / ih;

	cout << "real angle of ellpitc flow is                 : " << Psi[ih] << endl;
	cout << "Calculated Event angle of elliptic flow  is   : " << C_Psi[ih] << endl;
	// EP calculation done.


	TH1D *h_v2 = new TH1D("v2_histo","v2_histo",100, -1, 1);

	for( int i=0; i<N_Particle; i++){
		double dphi = phi[i] - Psi[ih] ; 
		h_v2->Fill(cos(ih*dphi));
	} // v2 calculation 

	v2_measured = h_v2->GetMean();


	cout << "real v2 values is 					: " << v2_input << endl;
    cout << "measured v2 value is				: " << v2_measured << endl;

}
